Machine learning to segment neutron images

Machine learning to segment neutron images

Anders Kaestner, Beamline scientist - Neutron Imaging

Laboratory for Neutron Scattering and Imaging
Paul Scherrer Institut

Lecture outline

In this lecture about machine learning to segment neutron images we will cover the following topics

1. Introduction
2. Limited data problem
3. Unsupervised segmentation
4. Supervised segmentation
5. Final problem: Segmenting root networks using convolutional NNs
6. Future Machine learning challenges in NI

Importing needed modules

This lecture needs some modules to run. We import all of them here.

import matplotlib.pyplot as plt
import seaborn as sn
import numpy as np
import pandas as pd
import skimage.filters as flt
import skimage.io as io
import matplotlib as mpl
from sklearn.cluster import KMeans
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import confusion_matrix
from matplotlib.colors import ListedColormap
from lecturesupport import plotsupport as ps
import pandas as pd
from sklearn.datasets import make_blobs
import scipy.stats as stats
import astropy.io.fits as fits
%matplotlib inline


from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg', 'png')
#plt.style.use('seaborn')
mpl.rcParams['figure.dpi'] = 300

import importlib
importlib.reload(ps);

Introduction

In this introduction part we give a starting point for the lecture. With topics like.

• Introduction to neutron imaging

  • Some words about the method
  • Contrasts

• Introduction to segmentation

  • What is segmentation
  • Noise and SNR

• Problematic segmentation tasks

  • Intro
  • Segmenation problems in neutron imaging

From these topics we will go into looking at how different machine learning techniques can be used to segment images. An in particular images obtained in neutron imaging experiments

What is an image?

In this lecture we are going to analyze the contents of images, but first lets define what an image is.

A very abstract definition:

• A pairing between spatial information (position)
• and some other kind of information (value).

In most cases this is a two- or three-dimensional position (x,y,z coordinates) and a numeric value (intensity)

In reality it can be many things like a picture, a radiograph, a matrix presentation of a raster scan and so on. In the case of volume images we are talking about volume representations like a tomography or time series like movies.

Science and Imaging

Images are great for qualitative analyses since our brains can quickly interpret them without large programming investements.

This is also the reason why image processing can be very frustrating. Seeing a feature in the image may be obvious to the human eye, but it can be very hard to write an algorithm that makes the same conclusion. In particular, when you are using traditional image processing techniques. This is an classic area for machine learning, these methods are designed and trained to recognize patterns and features in the images. Still, they can only perform well on the type of data they are trained with.

Proper processing and quantitative analysis is however much more difficult with images.

• If you measure a temperature, quantitative analysis is easy, $50K$.
• If you measure an image it is much more difficult and much more prone to mistakes, subtle setup variations, and confusing analyses

Furthermore in image processing there is a plethora of tools available

• Thousands of algorithms available
• Thousands of tools
• Many images require multi-step processing
• Experimenting is time-consuming

Some word about neutron imaging

Neutron imaging is a method based on transmission of the neutron radiation through a sample, i.e. the fundamental information is a radiograph. In this sense it is very much similar to the more known x-ray imaging. The basic outline in the figure below show the components source, collimator, sample, and detector. The neutron beam is slightly divergent, but it can mostly be approximated by a parallel beam.

```{figure} figures/collimator.pdf

scale: 75%

Schematic transmission imaging setup.

The intensity in a radiographic image proportional to the amount of radiation that remains after it was transmitted through the sample. The transmitted radiation is described by Beer-Lambert's law which in its basic form look like

$I=I_0\cdot{}e^{-\int_L \mu{}(x) dx}$
```{figure} figures/AttenuationLaw.pdf
---
scale: 30%
---
Transmission plot through a sample.

Where $\mu(x)$ is the attenuation coefficient of the sample at position x. This equation is a simplification as no source has a uniform spectrum and the attenuation coefficient depends on the radiation energy. We will touch the energy dispersive property in one of the segmentation examples later.

Single radiographs are relatively rare. In most experiments the radiographs are part of a time series aquisition or projections for the 3D tomography reconstuction. In this lecture we are not going very much into the details about the imaging method as such. This is a topic for other schools that are offered.

Neutron imaging contrast

The integral of Beer-Lambert's law turns into a sum for simple material compositions like in the sample in the figure below

```{figure} figures/porous_media_sand.pdf

scale: 50%

The transmission through a sample can be divied into contributions for each material.

The amount of radiation passing through the sample depends on the attenuation coefficients which are material and radiation specific. The periodic systems below illustrate how diffent elements attenuate for neutrons and X-rays

```{figure} figures/Periodic_table_neutron.pdf
---
scale: 75%
---
Attenuation coefficients for thermal neutrons.

```{figure} figures/Periodic_table_Xray.pdf

scale: 75%

Attenuation coefficients for X-rays. ```

Measurements are rarely perfect

There is no perfect measurement. This is also true for neutron imaging. The ideal image is the sample is distorted for many reasons. The figure below shows how an image of a sample can look after passing though the acquisition system. These quailty degradations will have an impact on the analysis of the image data. In some cases, it is possible to correct for some of these artifacts using classing image processing techniques. There are however also cases that require extra effort to correct the artifacts.

```{figure} figures/imperfect_imaging_system.pdf

scale: 75%

Schematic showing the error contributions in an imperfect imaging system. ```

Factors affecting the image quality

The list below provides some factors that affect the quality of the acquired images. Most of them can be handled by changing the imaging configuration in some sense. It can however be that the sample or process observed put limitiations on how much the acquisition can be tuned to obtain the perfect image.

• Resolution (Imaging system transfer functions)
• Noise
• Contrast
• Inhomogeneous contrast
• Artifacts

Resolution

The resolution is primarily determined optical transfer function of the imaging system. The actual resolution is given by the extents of the sample and how much the detector needs to capture in one image. This gives the field of view and given the number pixels in the used detector it is possible to calculate the pixel size. The pixel size limits the size of the smallest feature in the image that can be detected. The scintillator, which is used to convert neutrons into visible light, is chosen to

1. match the sampling rate given by the pixel size.
2. provide sufficient neutron capture to obtain sufficient light output for a given exposure time.

Noise

An imaging system has many noise sources, each with its own distribution e.g.

1. Neutron statistics - how many neutrons are collected in a pixel. This noise is Poisson distributed.
2. Photon statistics - how many photons are produced by each neutron. This noise is also Poisson distributed.
3. Thermal noise from the electronics which has a Gaussian distribution.
4. Digitation noise from converting the charges collected for the photons into digital numbers that can be transfered and stored by a computer, this noise has a binominal distribution.

The neutron statistics are mostly dominant in neutron imaging but in some cases it could also be that the photon statistics play a role.

Contrast

The contrast in the sample is a consequence of

1. how well the sample transmits the chosen radiation type. For neutrons you obtain good contrast from materials containing hydrogen or lithium while many metals are more transparent.
2. the amount of a specific element or material represented in a unit cell, e.g. a pixel (radiograph) or a voxel (tomography).

The objective of many experiments is to quantify the amount of a specific material. This could for example be the amount of water in a porous medium.

Good contrast between different image features is important if you want to segment them to make conclusions about the image content. Therefore, the radiation type should be chosen to provide the best contrast between the features.

Inhomogeneous contrast

The contrast in the raw radiograph depends much on the beam profile. These variations are easily corrected by normalizing the images by an open beam or flat field image.

• Biases introduced by scattering Scattering is the dominant interaction for many materials use in neutron imaging. This means that neutrons that are not passing straight though the sample are scattered and contribute to a background cloud of neutrons that build up a bias of neutron that are also detected and contribute to the

• Biases from beam hardening is a problem that is more present in x-ray imaging and is caused by that fact that the attenuation coefficient depends on the energy of the radiation. Higher energies have lower attenuation coefficient, thus will high energies penetrate the thicker samples than lower energies. This can be seen when a polychromatic beam is used.

Artifacts

Many images suffer from outliers caused by stray rays hitting the detector. Typical artefacts in tomography data are

• Lines, which are caused by outlier spots that only appear in single projections. These spot appear as lines in the reconstructed images.
• Rings are caused by stuck pixels which have the same value in a projections.

Introduction to segmentation

The basic task of segmentation is to identify regions of pixels with similar properties. This can be on different levels of abstraction. The lowest level is to create segments of pixels based on their pixel values like in the LANDSAT image below.

```{figure} figures/landsat_example.pdf

scale: 50%

Segmentation of a satelite image to identify different land types. ```

This type of segmentation is often done with the help of a histogram that shows the distribution of values in the image.

Different types of segmentation

• Semantic segmentation - pixel level

classifies all the pixels of an image into meaningful classes of objects. These classes are “semantically interpretable” and correspond to real-world categories. For instance, you could isolate all the pixels associated with a cat and color them green. This is also known as dense prediction because it predicts the meaning of each pixel.

• Instance segmentation - region level

Identifies each instance of each object in an image. It differs from semantic segmentation in that it doesn’t categorize every pixel. If there are three cars in an image, semantic segmentation classifies all the cars as one instance, while instance segmentation identifies each individual car.

Basic segmentation: Applying a threshold to an image

Start out with a simple image of a cross with added noise

$$ I(x,y) = f(x,y) $$

Here, we create a test image with two features embedded in uniform noise; a cross with values in the order of '1' and background with values in the order '0'. The figure below shows the image and its histogram. The histogram helps us to see how the graylevels are distributed which helps to decide where to put a threshold that segments the cross from the background.

fig,ax = plt.subplots(1,2,figsize=(7,3))
nx = 5; ny = 5;
xx, yy = np.meshgrid(np.arange(-nx, nx+1)/nx*2*np.pi, 
                     np.arange(-ny, ny+1)/ny*2*np.pi)
cross_im =   1.5*np.abs(np.cos(xx*yy))/(np.abs(xx*yy)+(3*np.pi/nx)) + np.random.uniform(-0.25, 0.25, size = xx.shape)
im=ax[0].imshow(cross_im, cmap = 'hot'); 
ax[1].hist(cross_im.ravel(),bins=10);

Applying a threshold to an image

By examining the image and probability distribution function, we can deduce that the underyling model is a whitish phase that makes up the cross and the darkish background

Applying the threshold is a deceptively simple operation

$$ I(x,y) = \begin{cases} 1, & f(x,y)\geq0.40 \\ 0, & f(x,y)<0.40 \end{cases}$$

threshold = 0.4; thresh_img = cross_im > threshold

fig,ax = plt.subplots(1,2,figsize=(8,4))
ax[0].imshow(cross_im, cmap = 'hot', extent = [xx.min(), xx.max(), yy.min(), yy.max()])
ax[0].plot(xx[np.where(thresh_img)]*0.9, yy[np.where(thresh_img)]*0.9,
           'ks', markerfacecolor = 'green', alpha = 0.5,label = 'Threshold', markersize = 15); ax[0].legend(fontsize=8);
ax[1].hist(cross_im.ravel(),bins=10); ax[1].axvline(x=threshold,color='r',label='Threshold'); ax[1].legend(fontsize=8);

Noise and SNR

Any measurement has a noise component and this noise has to be dealt with in some sense when the image is to be segmented. In the noise is not handled correctly it will cause many misclassified pixels. The noise strength of measured using the signal to noise ratio SNR. It is defined as the ratio between signal average and signal standard deviation.

The noise in neutron imaging mainly originates from the amount of captured neutrons. This noise is Poisson distributed and the signal to noise ratio is

$$SNR=\frac{E[x]}{s[x]}\sim\frac{N}{\sqrt{N}}=\sqrt{N}$$

where N is the number of captured neutrons. The figure below shows two neutron images acquired at 0.1s and 10s respectively. The plot shows the signal to noise ratio for different exposure times.

```{figure} figures/snrhanoi.pdf

scale: 75%

Signal to noise ratio for radiographies acqired with different exposure times. ```

Problematic segmentation tasks

Segmentation is rarely an obvious task. What you want to find is often obscured by other image features and unwanted artefacts from the experiment technique. If you take a glance at the painting by Bev Doolittle, you quickly spot the red fox in the middle. Looking a little closer at the painting, you'll recognize two indians on their spotted ponies. This example illustrates the problems you will encounter when you start to work with image segmentation.

```{figure} figures/doolittle_woodlandencounter.png

scale: 75%

Cases making the segmentation task harder than just applying a single thresshold. ```

Woodland Encounter Bev Doolittle

Typical image features that makes life harder

The basic segmentation shown in the example above can only be used under good conditions when the classes are well separated. Images from many experiments are unfortunately not well-behaved in many cases. The figure below shows four different cases when an advanced technique is needed to segment the image. Neutron images often show a combination of all cases.

```{figure} figures/trickycases.pdf

scale: 75%

Cases making the segmentation task harder than just applying a single thresshold. ```

The impact on the segmentation performance of all these cases can be reduced by proper experiment planning. Still, sometimes these improvements are not fully implemented in the experiment to be able to fulfill other experiment criteria.

Segmentation problems in neutron imaging

Limited data problem

Machine learning methods require a lot of training data to be able to build good models that are able to detect the features they are intended to.

Different types of limited data:

• Few data points or limited amounts of images

This is very often the case in neutron imaging. The number of images collected during an experiment session is often very low due to the long experiment duration and limited amount of beam time. This makes it hard to develop segmentation and analysis methods for single experiments. The few data points problem can partly be overcome by using data from previous experiment with similar characteristics.

• Unbalanced data

Unbalanced data means that the ratio between the data points with features you want detect and the total number data points is several orders of magnitude. E.g roots in a volume like the example we will look at later in this lecture. There is even a risk that the distribution of the wanted features is overlapped by the dominating background distribution.

```{figure} figures/classunbalance.pdf

scale: 100%

Two cases of unblanaced data; (a) the classes are well separated and the feature class is clearly visible in the tail distribution of the background and (b) the feature class is embeded in the background making it hard to detect. ``` Case (a) can most likely be segmented using one of the many histogram based thresholding methods proposed in literature.

• Little or missing training data

A complete set of training data contains both input data and labelled data. The input data is easy to obtain, it is the images you measured during your experiment. The labelled data is harder to get as it is a kind of chicken and egg problem. In particular, if your experiment data is limited. In that case, you would have to mark-up most of the available data to obtain the labeled data. Which doesn't make sense because

• then you'd already solved the task before even starting to train your segmentation algorithm.
• An algorithm based on learning doesn't improve the results, it only make it easier to handle large amounts of data.

Training data from NI is limited

In the introducing words about limited data essentially describes the problems that arise in neutron imaging.

• Long experiment times
• Few samples
• Some recycling from previous experiments is posible.

The experiment times are usually rather long which results in only few data sets per experiment. The advantage is that there are plenty experiments over the years that produce similar data which will add up to the pool of data.

Augmentation

Transfer learning

Unsupervised segmentation

Introducing clustering

With clustering methods you aim to group data points together into a limited number of clusters. Here, we start to look at an example where each data point has two values. The test data is generated using the make_blobs function.

test_pts = pd.DataFrame(make_blobs(n_samples=200, random_state=2018)[
                        0], columns=['x', 'y'])
plt.plot(test_pts.x, test_pts.y, 'r.');

The generated data set has two obvoius clusters, but if look closer it is even possible to identify three clusters. We will now use this data to try the k-means algorithm.

k-means

The k-means algorithm is one of the most used unsupervised clustering method. The user only have to provide the number of classes the algorithm shall find. Note: If you look for N groups you will always find N groups with K-Means, whether or not they make any sense

It is an iterative method that starts with a label image where each pixel has a random label assignment. Each iteration involves the following steps:

1. Compute current centroids based on the o current labels
2. Compute the value distance for each pixel to the each centroid. Select the class which is closest.
3. Reassign the class value in the label image.
4. Repeat until no pixels are updated.

The distance from pixel i to centroid j is usually computed as $||p_i - c_j||_2$.

It is important to note that k-means by definition is not position sensitive. The position can however be included as additional components in the data vectors.

k-means makes most sense to use on vector valued images where each pixel is represented by several values, e.g.:

1. Images from multimodal experiments like combined neutron and X-ray.
2. Wavelength resolved imaging

Basic clustering example

In this example we will use the blob data we previously generated and to see how k-means behave when we select different numbers of clusters.

fig, ax = plt.subplots(1,3,figsize=(15,4.5))

for i in range(3) :
    km = KMeans(n_clusters=i+2, random_state=2018); n_grp = km.fit_predict(test_pts)
    ax[i].scatter(test_pts.x, test_pts.y, c=n_grp)
    ax[i].set_title('{0} groups'.format(i+2))

When we select two clusters there is a natural separation between the two clusters we easily spotted by just looking at the data. When the number of clusters is increased to three, we again see a cluster separation that makes sense. Now, when the number of clusters is increased yet another time we see that one of the clusters is split once more. This time it is how ever questionable if the number of clusters makes sense. From this example, we see that it is important to be aware of problems related to over segmentation.

When can clustering be used on images?

Clustering applied to wavelength resolved imaging

The imaging techniques and its applications

The data

In this example we use a time of flight transmission spectrum from an assembly of six rectangular steel blocks produced by means of additive manufactoring. The typical approach to analyse this kind of information is to select some regions of interest and compute the average spectra from these regions. The average spectra are then analyzed using single Bragg edge fitting methods or Rietveld refinement.

Here, we will explore the posibility to identify regions with similar spectra using k-means clustering. The data is provided on the repository and has the dimensions 128x128x661, where the first to dimensions reflect the image size and the last is th number of time bins of the time of flight spectrum. This is a downsized version of the original data which has 512x512 image frames. Each frame has a relatively low signal to noise ratio, therefore we use the average image wtof to show the image contents.

tof  = np.load('../data/tofdata.npy')
wtof = tof.mean(axis=2)
plt.imshow(wtof);

Reshaping

The k-means requires vectors for each feature dimension, not images as we have in this data set. Therefore, we need to reshape our data. The reshaping is best done using the numpy array method reshape like this:

tofr=tof.reshape([tof.shape[0]*tof.shape[1],tof.shape[2]])
print("Input ToF dimensions",tof.shape)
print("Reshaped ToF data",tofr.shape)

Input ToF dimensions (128, 128, 661)
Reshaped ToF data (16384, 661)

The reshaped data tofr now has the dimensions 16384x661, i.e. each 128x128 pixel image has been replaced by a vector with the length 16384 elements. The number of time bins still remains the same. The reshapemethod is a cheap operation as it only changes the elements of the dimension vector and doesn't rearrange the data in memory.

Setting up and running k-means

When we set up k-means, we merely have to select the number of clusters. The choice depends on many factors. Particularly in this case with data points containing a large number of values we have a great degree of freedom in what could be considered the correct result. Thus we should ask ourselves what we expect from the clustering.

• We can clearly see that there is void on the sides of the specimens.
• There is also a separating band between the specimens.
• Finally we have to decide how many regions we want to find in the specimens. Let's start with two regions with different characteristics.

km = KMeans(n_clusters=4, random_state=2018)
c  = km.fit_predict(tofr).reshape(tof.shape[:2]) # Label image
kc = km.cluster_centers_.transpose()             # cluster centroid spectra

Results from the first try

fig,axes = plt.subplots(1,3,figsize=(18,5)); axes=axes.ravel()
axes[0].imshow(wtof,cmap='viridis'); axes[0].set_title('Average image')
p=axes[1].plot(kc);                  axes[1].set_title('Cluster centroid spectra'); axes[1].set_aspect(tof.shape[2], adjustable='box')
cmap=ps.buildCMap(p) # Create a color map with the same colors as the plot

im=axes[2].imshow(c,cmap=cmap); plt.colorbar(im);
axes[2].set_title('Cluster map');
plt.tight_layout()

This result was not very convincing...

We need more clusters

• Experiment data has variations on places we didn't expect k-means to detect as clusters.
• We need to increase the number of clusters!

What happens when we increase the unmber of clusters to ten?

km = KMeans(n_clusters=10, random_state=2018)
c  = km.fit_predict(tofr).reshape(tof.shape[:2]) # Label image
kc = km.cluster_centers_.transpose()             # cluster centroid spectra

Results of k-means with ten clusters

fig,axes = plt.subplots(1,3,figsize=(18,5)); axes=axes.ravel()
axes[0].imshow(wtof,cmap='viridis'); axes[0].set_title('Average image')
p=axes[1].plot(kc);                  axes[1].set_title('Cluster centroid spectra'); axes[1].set_aspect(tof.shape[2], adjustable='box')
cmap=ps.buildCMap(p) # Create a color map with the same colors as the plot

im=axes[2].imshow(c,cmap=cmap); plt.colorbar(im);
axes[2].set_title('Cluster map');
plt.tight_layout()

Interpreting the clusters

fig,axes = plt.subplots(1,2,figsize=(14,5)); axes=axes.ravel()
axes[0].matshow(np.corrcoef(kc.transpose()))
axes[1].plot(kc); axes[1].set_title('Cluster centroid spectra'); axes[1].set_aspect(tof.shape[2], adjustable='box')

Supervised segmentation

The supervised method needs to be told how to behave when it is presented with some new data. This is done during a training session.

```{figure} figures/teachingclass.pdf

scale: 100%

Supervised methods require a training session before performing the actual segmentation task. ```

The typical workflow with a supervised method involves the following three steps:

1. Training: Requires training data
2. Verification: Requires verification data
3. Inference: The images you want to segment

The data for the training and verification steps need the same type of data, i.e., data with markup that tell which parts of the image belongs to which segment type. It is how even not good if you use the same data for both training and verification. That would be a self-fulfilling prophecy and the performance numbers for the trained model will be misleadingly high. Typically, you want to have more data for training than for verification. After all, the training data is the part that builds the model and you want good statistic foundation for the model.

k nearest neighbors

The k nearest neighbors algorithm makes the inference based on a point cloud of training point. When the model is presented with a new point it computes the distance to the closest points in the model. The k in the algorithm name indicates how many neighbors should be considered. E.g. k=3 means that the major class of the three nearest neighbors is assigned the tested point. Looking at the example below we would say that using three neighbors the

• Green hiker would claim he is in a spruce forest.
• Orange hiker would claim he is in a mixed forest.
• Blue hiker would claim he is in a birch forest.

```{figure} figures/forest.pdf

scale: 100%

Depending on the where the hiker is standing he makes the the conclusion that is either in a birch or spruce forest. ```

Create example data for supervised segmentation

blob_data, blob_labels = make_blobs(n_samples=100, random_state=2018)
test_pts = pd.DataFrame(blob_data, columns=['x', 'y'])
test_pts['group_id'] = blob_labels
plt.scatter(test_pts.x, test_pts.y, c=test_pts.group_id, cmap='viridis');

Detecting unwanted outliers in neutron images

orig= fits.getdata('../data/spots/mixture12_00001.fits')
annotated=io.imread('../data/spots/mixture12_00001.png'); mask=(annotated[:,:,1]==0)
r=600; c=600; w=256
ps.magnifyRegion(orig,[r,c,r+w,c+w],[15,7],vmin=400,vmax=4000,title='Neutron radiography')

Marked-up spots

```{figure} figures/markedspots.pdf

scale: 100%

Two cases of unblanaced data; (a) the classes are well separated and the feature class is clearly visible in the tail distribution of the background and (b) the feature class is embeded in the background making it hard to detect. ```

Baseline - Traditional spot cleaning algorithm

```{figure} figures/spotclean_algorithm.pdf

scale: 100%

Schematic description of the spot cleaning baseline algorithm ```

Parameters

• N Width of median filter.
• k Threshold level for outlier detection.

The spot cleaning algorithm

The baseline algorithm is here implemented as a python function that we will use when we compare the performance of the CNN algorithm. This is the most trivial algorithm for spot cleaning and there are plenty other algorithms to solve this task.

def spotCleaner(img, threshold=0.95, selem=np.ones([3,3])) :
    fimg=img.astype('float32')
    mimg = flt.median(fimg,selem=selem)
    timg = threshold < np.abs(fimg-mimg)
    cleaned = mimg * timg + fimg * (1-timg)
    return (cleaned,timg)

baseclean,timg = spotCleaner(orig,threshold=1000)
ps.magnifyRegion(baseclean,[r,c,r+w,c+w],[12,3],vmin=400,vmax=4000,title='Cleaned image')
ps.magnifyRegion(timg,[r,c,r+w,c+w],[12,3],vmin=0,vmax=1,title='Detection image')

k nearest neighbors to detect spots

selem=np.ones([3,3])
forig=orig.astype('float32')
mimg = flt.median(forig,selem=selem)
d = np.abs(forig-mimg)

fig,ax=plt.subplots(1,1,figsize=(8,5))
h,x,y,u=ax.hist2d(forig[:1024,:].ravel(),d[:1024,:].ravel(), bins=100);
ax.imshow(np.log(h[::-1]+1),vmin=0,vmax=3,extent=[x.min(),x.max(),y.min(),y.max()])
ax.set_xlabel('Input image - $f$'),ax.set_ylabel('$|f-med_{3x3}(f)|$'),ax.set_title('Log bivariate histogram');

Prepare data

Training data

trainorig = forig[:,:1000].ravel()
traind    = d[:,:1000].ravel()
trainmask = mask[:,:1000].ravel()

train_pts = pd.DataFrame({'orig': trainorig, 'd': traind, 'mask':trainmask})

Test data

testorig = forig[:,1000:].ravel()
testd    = d[:,1000:].ravel()
testmask = mask[:,1000:].ravel()

test_pts = pd.DataFrame({'orig': testorig, 'd': testd, 'mask':testmask})

Train the model

k_class = KNeighborsClassifier(1)
k_class.fit(train_pts[['orig', 'd']], train_pts['mask']) 

KNeighborsClassifier(n_neighbors=1)

Inspect decision space

xx, yy = np.meshgrid(np.linspace(test_pts.orig.min(), test_pts.orig.max(), 100),
                     np.linspace(test_pts.d.min(), test_pts.d.max(), 100),indexing='ij');
grid_pts = pd.DataFrame(dict(x=xx.ravel(), y=yy.ravel()))
grid_pts['predicted_id'] = k_class.predict(grid_pts[['x', 'y']])
plt.scatter(grid_pts.x, grid_pts.y, c=grid_pts.predicted_id, cmap='gray'); plt.title('Testing Points'); plt.axis('square');

Apply knn to unseen data

pred = k_class.predict(test_pts[['orig', 'd']])
pimg = pred.reshape(d[1000:,:].shape)

fig,ax = plt.subplots(1,3,figsize=(15,6))
ax[0].imshow(forig[1000:,:],vmin=0,vmax=4000), ax[0].set_title('Original image')
ax[1].imshow(pimg), ax[1].set_title('Predicted spot')
ax[2].imshow(mask[1000:,:]),ax[2].set_title('Annotated spots');

Performance check

cmbase = confusion_matrix(mask[:,1000:].ravel(), timg[:,1000:].ravel(), normalize='all')
cmknn  = confusion_matrix(mask[:,1000:].ravel(), pimg.ravel(), normalize='all')

fig,ax = plt.subplots(1,2,figsize=(10,4))
sn.heatmap(cmbase, annot=True,ax=ax[0]), ax[0].set_title('Confusion matrix baseline')
sn.heatmap(cmknn, annot=True,ax=ax[1]), ax[1].set_title('Confusion matrix k-NN')

(<AxesSubplot:title={'center':'Confusion matrix k-NN'}>,
 Text(0.5, 1.0, 'Confusion matrix k-NN'))

Some remarks about k-nn

• It takes more time to process
• You need to prepare training data
  • Annotation takes time...
  • Here we used the segmentation on the same type of image
  • We should normalize the data
  • This was a raw projection, what happens if we use a flat field corrected image?
• Finds more spots than baseline
• Data is very unbalanced, try a selection of non-spot data for training.
  • Is it faster?
  • Is there a drop segmentation performance?

Note There are other spot detection methods that perform better than the baseline.

Convolutional neural networks for segmentation

Training data

We have two choices:

1. Use real data
   • requires time consuming markup to provide training data
   • corresponds to real life images
2. Synthesize data
   • flexible and provides both 'dirty' data and ground truth.
   • model may not behave as real data

Preparing real data

We will use the spotty image as training data for this example

Split image to tiles

The images we get are usually much larger than is feasible to handle for a CNN with given hardware. Therefore we will split the image into many smaller images - tiles. The function below will do this job for us and place the tiles in a list.

```{figure} figures/WorkflowWithValidationSet.pdf

scale: 100%

How the three data sets training, validation, and test are used when a network is trained and optimized. ```

def splitTiles(img,size=[64,64]) :
    dims = img.shape
    nTiles = [dims[0]//(size[0]), dims[0]//(size[0])]
    
    tiles = []
    
    for x in range(nTiles[0]) :
        for y in range(nTiles[1]) :
            tiles.append(img[x*size[0]:(x+1)*size[0],y*size[1]:(y+1)*size[1]])
            
    return tiles

Now we apply the tile splitter to our original image and the mask image where the spots are annotated.

origTiles = splitTiles(orig);
maskTiles = splitTiles(mask);

Lets inspect some tiles

fig,ax = plt.subplots(2,5,figsize=(15,8))
ax=ax.ravel()
for idx,item in enumerate(np.random.randint(len(origTiles), size=5)) :
    ax[idx].imshow(origTiles[item],vmin=200,vmax=4000)
    ax[idx+5].imshow(maskTiles[item])

Prepare training, validation, and test data

Any analysis system must be verified to be demonstrate its performance and to further optimize it.

For this we need to split our data into three categories:

1. Training data
2. Test data
3. Validation data

The tile splitting produced a list of small images. Now we need to divide the list into three parts before we can start training our network.

```{figure} figures/WorkflowWithValidationSet.pdf

scale: 100%

How the three data sets training, validation, and test are used when a network is trained and optimized. ```

Training	Validation	Test
70%	15%	15%

Build a CNN for spot detection and cleaning

We need:

• Data
• Tensorflow
  • Data provider
  • Model design

#### Data provider

Tensorflow needs a

#### Model design

Segmenting root networks in the rhizosphere using an U-Net

Background

• Soil and in particular the rhizosphere are of central interest for neutron imaging users.
• The experiments aim to follow the water distribution near the roots.
• The roots must be identified in 2D and 3D data
• Today: much of this mark-up is done manually!

Available data

Considered NN models

Loss functions

Training

Results

Summary

Future Machine learning challenges in neutron imaging

Concluding remarks